plot_ly() and ggplotly() functionsplot_geo()We will work with COVID data downloaded from the New York Times. The dataset consists of COVID-19 cases and deaths in each US state during the course of the COVID epidemic.
The objective of this lab is to explore relationships between cases, deaths, and population sizes of US states, and plot data to demonstrate this
# Combine data by state
cv_states_readin <- as.data.frame(data.table::fread("https://raw.githubusercontent.com/nytimes/covid-19-data/master/us-states.csv") )
state_pops <- as.data.frame(data.table::fread("https://raw.githubusercontent.com/COVID19Tracking/associated-data/master/us_census_data/us_census_2018_population_estimates_states.csv"))
state_pops$abb <- state_pops$state
state_pops$state <- state_pops$state_name
state_pops$state_name <- NULL
cv_states <- merge(cv_states_readin, state_pops, by="state")
head, and tail of the datadim(cv_states)
## [1] 39114 9
head(cv_states)
## state date fips cases deaths geo_id population pop_density abb
## 1 Alabama 2021-08-21 1 659750 12000 1 4887871 96.50939 AL
## 2 Alabama 2021-09-17 1 764839 13048 1 4887871 96.50939 AL
## 3 Alabama 2020-07-24 1 76005 1438 1 4887871 96.50939 AL
## 4 Alabama 2021-09-03 1 714860 12394 1 4887871 96.50939 AL
## 5 Alabama 2020-03-18 1 51 0 1 4887871 96.50939 AL
## 6 Alabama 2022-02-09 1 1254032 17452 1 4887871 96.50939 AL
tail(cv_states)
## state date fips cases deaths geo_id population pop_density abb
## 39109 Wyoming 2021-04-09 56 56873 701 56 577737 5.950611 WY
## 39110 Wyoming 2021-09-18 56 83958 918 56 577737 5.950611 WY
## 39111 Wyoming 2020-06-01 56 910 17 56 577737 5.950611 WY
## 39112 Wyoming 2021-04-20 56 57456 705 56 577737 5.950611 WY
## 39113 Wyoming 2020-04-11 56 343 0 56 577737 5.950611 WY
## 39114 Wyoming 2021-01-06 56 45890 464 56 577737 5.950611 WY
str(cv_states)
## 'data.frame': 39114 obs. of 9 variables:
## $ state : chr "Alabama" "Alabama" "Alabama" "Alabama" ...
## $ date : IDate, format: "2021-08-21" "2021-09-17" ...
## $ fips : int 1 1 1 1 1 1 1 1 1 1 ...
## $ cases : int 659750 764839 76005 714860 51 1254032 159169 932250 3953 104786 ...
## $ deaths : int 12000 13048 1438 12394 0 17452 2558 16503 114 1882 ...
## $ geo_id : int 1 1 1 1 1 1 1 1 1 1 ...
## $ population : int 4887871 4887871 4887871 4887871 4887871 4887871 4887871 4887871 4887871 4887871 ...
## $ pop_density: num 96.5 96.5 96.5 96.5 96.5 ...
## $ abb : chr "AL" "AL" "AL" "AL" ...
cv_states$date <- as.Date(cv_states$date, format="%Y-%m-%d")
state_list <- unique(cv_states$state)
cv_states$state <- factor(cv_states$state, levels = state_list)
abb_list <- unique(cv_states$abb)
cv_states$abb <- factor(cv_states$abb, levels = abb_list)
cv_states = cv_states[order(cv_states$state, cv_states$date),]
str(cv_states)
## 'data.frame': 39114 obs. of 9 variables:
## $ state : Factor w/ 52 levels "Alabama","Alaska",..: 1 1 1 1 1 1 1 1 1 1 ...
## $ date : Date, format: "2020-03-13" "2020-03-14" ...
## $ fips : int 1 1 1 1 1 1 1 1 1 1 ...
## $ cases : int 6 12 23 29 39 51 78 106 131 157 ...
## $ deaths : int 0 0 0 0 0 0 0 0 0 0 ...
## $ geo_id : int 1 1 1 1 1 1 1 1 1 1 ...
## $ population : int 4887871 4887871 4887871 4887871 4887871 4887871 4887871 4887871 4887871 4887871 ...
## $ pop_density: num 96.5 96.5 96.5 96.5 96.5 ...
## $ abb : Factor w/ 52 levels "AL","AK","AZ",..: 1 1 1 1 1 1 1 1 1 1 ...
head(cv_states)
## state date fips cases deaths geo_id population pop_density abb
## 571 Alabama 2020-03-13 1 6 0 1 4887871 96.50939 AL
## 380 Alabama 2020-03-14 1 12 0 1 4887871 96.50939 AL
## 431 Alabama 2020-03-15 1 23 0 1 4887871 96.50939 AL
## 33 Alabama 2020-03-16 1 29 0 1 4887871 96.50939 AL
## 168 Alabama 2020-03-17 1 39 0 1 4887871 96.50939 AL
## 5 Alabama 2020-03-18 1 51 0 1 4887871 96.50939 AL
tail(cv_states)
## state date fips cases deaths geo_id population pop_density abb
## 38547 Wyoming 2022-03-18 56 155907 1769 56 577737 5.950611 WY
## 39053 Wyoming 2022-03-19 56 155907 1769 56 577737 5.950611 WY
## 38906 Wyoming 2022-03-20 56 155907 1769 56 577737 5.950611 WY
## 38959 Wyoming 2022-03-21 56 155907 1769 56 577737 5.950611 WY
## 39020 Wyoming 2022-03-22 56 155988 1783 56 577737 5.950611 WY
## 38587 Wyoming 2022-03-23 56 155988 1783 56 577737 5.950611 WY
head(cv_states)
## state date fips cases deaths geo_id population pop_density abb
## 571 Alabama 2020-03-13 1 6 0 1 4887871 96.50939 AL
## 380 Alabama 2020-03-14 1 12 0 1 4887871 96.50939 AL
## 431 Alabama 2020-03-15 1 23 0 1 4887871 96.50939 AL
## 33 Alabama 2020-03-16 1 29 0 1 4887871 96.50939 AL
## 168 Alabama 2020-03-17 1 39 0 1 4887871 96.50939 AL
## 5 Alabama 2020-03-18 1 51 0 1 4887871 96.50939 AL
summary(cv_states)
## state date fips cases
## Washington : 793 Min. :2020-01-21 Min. : 1.00 Min. : 1
## Illinois : 790 1st Qu.:2020-09-05 1st Qu.:16.00 1st Qu.: 49810
## California : 789 Median :2021-03-12 Median :29.00 Median : 209010
## Arizona : 788 Mean :2021-03-12 Mean :29.78 Mean : 536429
## Massachusetts: 782 3rd Qu.:2021-09-16 3rd Qu.:44.00 3rd Qu.: 665336
## Wisconsin : 778 Max. :2022-03-23 Max. :72.00 Max. :9079164
## (Other) :34394
## deaths geo_id population pop_density
## Min. : 0.0 Min. : 1.00 Min. : 577737 Min. : 1.292
## 1st Qu.: 884.2 1st Qu.:16.00 1st Qu.: 1805832 1st Qu.: 43.659
## Median : 3525.5 Median :29.00 Median : 4468402 Median : 107.860
## Mean : 8812.6 Mean :29.78 Mean : 6419723 Mean : 422.717
## 3rd Qu.:10688.8 3rd Qu.:44.00 3rd Qu.: 7535591 3rd Qu.: 229.511
## Max. :88461.0 Max. :72.00 Max. :39557045 Max. :11490.120
## NA's :741
## abb
## WA : 793
## IL : 790
## CA : 789
## AZ : 788
## MA : 782
## WI : 778
## (Other):34394
min(cv_states$date)
## [1] "2020-01-21"
max(cv_states$date)
## [1] "2022-03-23"
new_cases and new_deaths and correct outliersnew_cases, and new deaths, new_deaths:
new_cases equal to the difference between cases on date i and date i-1, starting on date i=2for (i in 1:length(state_list)) {
cv_subset = subset(cv_states, state == state_list[i])
cv_subset = cv_subset[order(cv_subset$date),]
# add starting level for new cases and deaths
cv_subset$new_cases = cv_subset$cases[1]
cv_subset$new_deaths = cv_subset$deaths[1]
for (j in 2:nrow(cv_subset)) {
cv_subset$new_cases[j] = cv_subset$cases[j] - cv_subset$cases[j-1]
cv_subset$new_deaths[j] = cv_subset$deaths[j] - cv_subset$deaths[j-1]
}
cv_states$new_cases[cv_states$state==state_list[i]] = cv_subset$new_cases
cv_states$new_deaths[cv_states$state==state_list[i]] = cv_subset$new_deaths
}
# Focus on recent dates
cv_states <- cv_states %>% dplyr::filter(date >= "2021-07-01")
ggplotly for EDA: See if there are outliers or values that don’t make sense for new_cases and new_deaths. Which states and which dates have strange values?p1 <- ggplot(cv_states,aes(x=date,y=new_cases, colour=state)) +
geom_line() +
ylab("New Cases") +
xlab("Date") +
labs(title="New COVID Cases Overtime by State") +
theme_minimal()
ggplotly(p1)
p2 <- ggplot(cv_states, aes(x=date,y=new_deaths,colour=state)) +
geom_line() +
ylab("New Deaths") +
xlab("Date") +
labs(title="New COVID Deaths Overtime by State") +
theme_minimal()
ggplotly(p2)
new_cases or new_deaths to 0It was strange that there are negative deaths and cases, which were results of redefining terms by state authorities.
cv_states$new_cases[cv_states$new_cases < 0] = 0
cv_states$new_deaths[cv_states$new_deaths < 0] = 0
p3 <- ggplot(cv_states, aes(x=date, y=new_deaths, color = state)) +
geom_line() +
geom_point(size=0.5,alpha=0.5) +
ylab("New Deaths") +
xlab("Date") +
labs(title="New COVID Deaths Overtime by State (outliers removed)")+
theme_minimal()
ggplotly(p3)
Add population-normalized (by 100,000) variables for each variable type (rounded to 1 decimal place). Make sure the variables you calculate are in the correct format (numeric). You can use the following variable names:
per100k = cases per 100,000 populationnewper100k= new cases per 100,000deathsper100k = deaths per 100,000newdeathsper100k = new deaths per 100,000Add a “naive CFR” variable representing deaths / cases on each date for each state
Create a dataframe representing values on the most recent date, cv_states_today
cv_states$per100k = as.numeric(format(round(cv_states$cases/(cv_states$population/100000),1),nsmall=1))
cv_states$newper100k = as.numeric(format(round(cv_states$new_cases/(cv_states$population/100000),1),nsmall=1))
cv_states$deathsper100k = as.numeric(format(round(cv_states$deaths/(cv_states$population/100000),1),nsmall=1))
cv_states$newdeathsper100k = as.numeric(format(round(cv_states$new_deaths/(cv_states$population/100000),1),nsmall=1))
cv_states = cv_states %>% mutate(naive_CFR = round((deaths*100/cases),2))
max_date <- max(cv_states$date)
cv_states_today = cv_states %>% filter(date==as.Date(max_date))
plot_ly()plot_ly() representing pop_density vs. various variables (e.g. cases, per100k, deaths, deathsper100k) for each state on most recent date (cv_states_today)
# Restrict sizes to a range, specify how we want the markers to be
cv_states_today %>%
plot_ly(x=~pop_density, y=~cases,type="scatter", mode="markers", color=~state,
size=~population, sizes=c(5,70),marker=list(sizemode='diameter'),
opacity=0.5) %>%
layout(title = 'Cases by Population Density and State', xaxis = list(title = 'Population Density'),
yaxis = list(title = 'Cases'), legend = list(title=list(text='<b> State </b>')))
# Remove the outlier not because it's an error, but just to make observing the other states' data easier
cv_states_today %>%
filter(state != "District of Columbia") %>%
plot_ly(x=~pop_density, y=~cases, color=~state, type="scatter", mode="markers",
size=~population, sizes=c(5,70), marker=list(sizemode="diameter", opacity=0.5)) %>%
layout(title = 'Cases by Population Density and State', xaxis = list(title = 'Population Density'),
yaxis = list(title = 'Cases'), legend = list(title=list(text='<b> State </b>')))
hovermode = "compare"# Paste0 puts spaces in between the categories
cv_states_today %>%
filter(state != "District of Columbia") %>%
plot_ly(x=~pop_density, y=~cases, color=~state, type="scatter", mode="markers",
size=~population, sizes=c(5,70),
marker=list(sizemode="diameter", opacity=0.5), hover_info="text",
text=~paste(
paste0("State: ", state),
paste0("Cases per 100k: ", per100k),
paste("Death per 100k:", deathsper100k),
sep="<br>")) %>%
layout(title = 'Population-normalized Cases per 100k',
xaxis = list(title = 'Population Density'),
yaxis = list(title = 'Cases per 100k'),
hovermode="compare",
legend = list(title=list(text='<b> State </b>')))
ggplotly()pop_density vs. newdeathsper100k create a chart with the same variables using ggplotly()pop_density correlates with newdeathsper100k?p4 <- cv_states_today %>%
filter(state != "District of Columbia") %>%
ggplot(aes(x=pop_density, y=newdeathsper100k, colour = state, size=population)) +
geom_point() +
ylab("New Deaths per 100k") +
xlab("Population Density") +
labs(title="Population Density versus New Deaths per 100k") +
theme_minimal()
ggplotly(p4)
naive_CFR for all states over time using plot_ly()
naive_CFR for the states that had an increase in September. How have they changed over time?new_cases and new_deaths together in one plot. Hint: use add_layer()
cv_states %>%
plot_ly(x=~date, y=~naive_CFR, color=~state, type="scatter", mode="lines") %>%
layout(title = 'Naive CFR Overtime',
xaxis = list(title = 'Date'),
yaxis = list(title = 'Naive CFR'),
legend = list(title=list(text='<b> State </b>')))
cv_states %>%
filter(state=="Florida") %>%
plot_ly(x=~date, y=~new_cases, type="scatter", mode = "lines",
name = "Cases") %>%
add_lines(x=~date, y=~new_deaths, type="scatter", mode="lines",
name="Deaths") %>%
layout(title = 'Florida\'s COVID Cases and Deaths Overtime',
xaxis = list(title = 'Date'),
yaxis = list(title = 'New Cases'),
legend = list(title=list(text='<b> Type </b>')))
Create a heatmap to visualize new_cases for each state on each date greater than July 1st, 2021 - Start by mapping selected features in the dataframe into a matrix using the tidyr package function pivot_wider(), naming the rows and columns, as done in the lecture notes - Use plot_ly() to create a heatmap out of this matrix. Which states stand out?
cv_states_mat <- cv_states %>%
select(state, date, new_cases) %>%
filter(date > "2021-07-01")
# Change pivot so that we have date as a column
cv_states_mat2 <- as.data.frame(pivot_wider(cv_states_mat, names_from=state,
values_from=new_cases))
# Here, colnames is not a variable so we don't need to write a function
cv_states_mat2 <- cv_states_mat2 %>%
column_to_rownames("date") %>%
as.matrix()
# Z value for color intensity
plot_ly(x=colnames(cv_states_mat2), y=rownames(cv_states_mat2), z=~cv_states_mat2, type="heatmap") %>%
layout(title = 'Heatmap for New Cases by State Overtime',
xaxis = list(title = 'State'),
yaxis = list(title = 'Date'))
Lab 10b questions 1-2, lab 11 questions 0-10. Upload html or pdf for both lab Rmd’s to quercus.